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Averaging and evolving inhomogeneities are non-commuting operations. This implies the exis- 
tence of deviations of an averaged model from the standard Friedmann-Lemaitre cosmologies. We 
quantify these deviations, encoded in a backreaction parameter, in the framework of Newtonian 
cosmology. We employ the linear theory of gravitational instability in the Eulerian and Lagrangian 
approaches, as well as the spherically- and plane-symmetric solutions as standards of reference. We 
propose a model for the evolution of the average characteristics of a spatial domain for generic initial 
conditions that contains the spherical top-hat model and the planar collapse model as exact sub 
cases. A central result is that the backreaction term itself, calculated on sufficiently large domains, 
is small but, still, its presence can drive the cosmological parameters on the averaging domain far 
away from their global values of the standard model. We quantify the variations of these parameters 
in terms of the fluctuations in the initial data as derived from the power spectrum of initial cold 
dark matter density fluctuations. E.g. in a domain with a radius of lOOMpc today and initially 
one-cr fluctuations, the density parameters deviate from their homogeneous values by 15%; three-cr 
fluctuations lead to deviations larger than 100%. 
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I. INTRODUCTION 

The standard picture of the evolution of the Universe 
relates typical observables to global scalar cosmological 
parameters such as the mean matter density and the 
global expansion rate. These parameters are derived 
from solutions of Einstein's or Newton's laws of gravita- 
tionally interacting homogeneous systems. Moreover, the 
class of homogeneous systems is usually narrowed to the 
isotropic ones, in most cases to the simplest, spatially flat 
model, first discussed by Einstein and de Sitter at the be- 
ginning of this century. This standard model has survived 
in spite of the tremendous development of our knowledge 
about inhomogeneities in the Universe. The key argu- 
ment of applicability of the homogeneous-isotropic mod- 
els is the conjecture that on average the Universe may 
still follow the evolution of these simple models. Without 
specifying the notion of averaging and without entering a 
sophisticated discussion of how one should average an in- 
homogeneous spacetime, it can be definitely said that this 
conjecture is in itself courageous given the nonlinearity 
of Einstein's or Newton's laws, the long-range attractive 
forces, and the fact that the solutions are linearly unsta- 
ble to inhomogeneous perturbations in the matter vari- 
ables. Because of the continued success of the standard 
model as an explanation of a variety of orthogonal obser- 
vations, the need to replace the homogeneous-isotropic 
models by averaged inhomogeneous models is not obvious 
to most researchers in cosmology (from a physical point 
of view it should be obvious). Notwithstanding we will 
add a substantial contribution to supporting this need in 
the present article, shedding new light on the interpre- 
tation of observations and N-body simulations in spatial 
domains that cover only a few percent of the Hubble vol- 
ume. 

Recent progress in the understanding of the effective 
(i.e. spatially averaged) dynamical properties of inho- 
mogeneous cosmological models may be summarized by 
saying that one can unambiguously replace the stan- 
dard Friedmann equations governing the homogeneous- 
isotropic models by corresponding equations for the spa- 
tially averaged variables, provided we focus our attention 
to scalar dynamical variables. The averaged equations 
feature an additional source, we call it "backreaction" 
of the inhomogeneities on the average expansion, which 
will be the basis of our investigations. The restriction to 
scalars has its root in the averaging problem of general 
relativity jlj , where a straightforward answer to the ques- 
tion of how one should average an inhomogeneous space- 
time metric is not at hand. For scalars spatial averaging 
can be properly defined, provided we have a foliation 
of spacetime. The results of averaging the Newtonian 
equations for the evolution of a single dust component 
P], where also the choice of foliation is not a problem, 
form the basis of the current work. The corresponding 
equations in general relativity || will also be put into 
perspective. 



We shall not review or relate other numerous work on 
the "backreaction problem", partly because we want to 
concentrate on a quantitative investigation of the general- 
ized Friedmann equations , and partly because the no- 
tion of "backreaction" as we employ it differs from other 
notions; for example, many authors study the averaging 
problem perturbatively in general relativity and consider 
as "backreaction" any deviation in curvature from the flat 
Friedmann cosmologies. The approximations employed 
sometimes even imply that the perturbed models remain 
within the class of Friedmann-Lemaitre cosmologies. To 
review the different averaging procedures was recently at- 
tempted by Stoeger et al. to review the vastly differ- 
ing methods of how the "backreaction" can be estimated 
is a considerable, but desirable program which is also 
complicated by different gauge choices in a relativistic 
treatment. We will mention some of these works in the 
main text mainly to emphasize that the "backreaction" 
we are talking about is not a genuinely relativistic effect, 
but is naturally present in the Newtonian framework as 
well. 

Also, in the present work, we shift the emphasis to 
fluctuations of the matter variables as a result of backre- 
action, and we carry out this investigation within a stan- 
dard cosmological model in which there is no global con- 
tribution to backreaction. In earlier, general relativistic 
investigations researchers mainly attempted to calculate 
this global effect. We shall quantify the influence of the 
backreaction on the expansion in finite domains, small 
compared to the horizon, using the Newtonian approxi- 
mation. Consistently, we have to assume that on some 
(very large) scale the Universe expands according to the 
Friedmann equation. Possible boundary conditions that 
respect the cosmological principle are rare and practically 
restricted to be periodic in which case the boundary is 
empty and thus the backreaction vanishes globally on the 
periodicity scale (compare J|], especially Appendix A). 
We will assume that this scale is of the order of the Hub- 
ble radius. Hence, with our Newtonian treatment, we are 
not able to give any quantitative global results, but we 
shall see that even domains with a size of hundreds of 
Mpc's are influenced by backreaction. 

Both N-body simulations and most analytical approxi- 
mations using perturbation theory rely on the Newtonian 
approximation to describe the formation of structure in 
the Universe. The initial conditions are often specified 
using a Gaussian random field for the density contrast. 
These two assumptions also enter our calculations of the 
backreaction. Note that both the numerical and analyti- 
cal approaches enforce a globally vanishing backreaction 
by imposing periodic boundary conditions on the inho- 
mogeneities. However, contrary to current N-body sim- 
ulations, we are not limited by box size and resolution, 
but we are limited by the validity of Zel'dovich's approx- 
imation, which we employ as evolution model. The latter 
has been tested thoroughly in comparison with N-body 
simulations and we shall add another test concerning its 
average performance. 
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We view the present work as a dynamical approach 
to cosmic variance: not only the density contrast, but 
also fluctuations in the peculiar-velocity gradient (e.g. 
fluctuations in shear, expansion and vorticity) are quan- 
tified by their effective impact on arbitrary Newtonian 
portions of the Universe. On such portions the inter- 
play of backreaction (condensed into an additional "cos- 
mological parameter") with the standard parameters of 
the "cosmic triangle" (matter density parameter, "curva- 
ture" parameter and cosmological constant, ||) implies 
scale-dependence of the fluctuations of cosmological pa- 
rameters. 

This work may be approached with the following guide- 
lines: 

• The "cosmic triangle" fluctuates on a given spa- 
tial scale. Comparing observations on a given 
scale with cosmological parameters of the standard 
model has to be subjected to statistical uncertain- 
ties of the parameters due to backreaction. 

• The influence of backreaction may be two-fold: 
for dominating shear fluctuations the effective den- 
sity causing the expansion is larger than the ac- 
tual matter density, thus mimicking a "kinemati- 
cal dark matter" ; for dominating expansion fluctu- 
ations backreaction acts accelerating, thus mimick- 
ing the effect of a (positive) cosmological constant. 

• Spatial portions of the Universe may collapse as a 
result of backreaction, even without the presence of 
over-densities. Averaged generic inhomogeneities 
invoke deviations in the time-dependence of the 
scale factor from spherical or plane collapse. 

• Although in some cases the backreaction source 
term itself may be numerically small, its impact 
can be large, since already a small source term is 
capable of driving the average system far away from 
the homogeneous solutions. 

We have chosen to organize the investigation of the ef- 
fect in the following way. In Sect. [n]we briefly review the 
generalized Friedmann equations following from averag- 
ing inhomogeneous dust cosmologies in the Newtonian 
framework and provide different representations of the 
backreaction that are used in this article. Wc then move 



in Sect. Ill to dynamical models that allow estimating the 
backreaction term: the Eulerian linear perturbation the- 
ory, and the Lagrangian linear perturbation theory (re- 
stricted to Zel'dovich's approximation) serve as approx- 
imate models. The plane and spherical collapse models 
serve as exact reference solutions. 



We then illustrate the results in Sect. IV in terms of the 
time evolution of cosmological parameters: we examine 
the evolution of the scale factor for typical initial condi- 
tions; then we explore derived parameters: the Bubble- 
parameter and the deceleration parameter as well as the 
cosmological density parameters. In Sect. [v| we discuss 



some related issues: we put the general relativistic equa- 
tions into perspective, and we explain the role of current 
N-body simulations in view of our results. 

Technicalities are left to appendices: Appendix ^ gives 
a short derivation of the generalized expansion law, Ap- 
pendix ^ lists useful writings of the principal scalar in- 
variants of a tensor field, and in Appendix |^ we relate 
the fluctuations of spatial invariants, entering the "back- 
reaction' term, to the initial power spectrum of the inho- 
mogeneities in the density contrast. 



II. AVERAGED DYNAMICS IN THE 
NEWTONIAN APPROXIMATION 

In this section we recall results obtained previously by 
Buchert and Ehlers |2j and summarize alternative writ- 
ings of the backreaction source term. 



A. The generalized Friedmann equation 

In the Newtonian approximation the expansion of a 
domain is influenced by the inhomogeneities inside the 
domain S. 

We shall investigate fields averaged over a simply- 
connected spatial domain T> at time t which evolved out 
of the initial domain T>i at time t\, conserving the mass 
inside. The averaged scale factor a-p, depending on con- 
tent, shape and position of the spatial domain of aver- 
aging T>, is defined via the domain's volume V(t) = \T>\ 
and the initial volume V\ = V(ti) = \T>i\: 



(1) 



In Fig. |l| we illustrate the evolution of an initial domain 
T>i with the Bubble flow, resulting in a simple rescaling 
by the global scale factor a(t). In a more general inho- 
mogeneous setting the scale factor a-p(t) is defined via 
Eq. (Eh. 



a(t)Oi 




D(t) 




FIG. 1. The evolution of an initial domain T>\ within 
the Hubble flow with global scale factor a(t) (left) and for a 
general inhomogeneous setting (right). 
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With (-) v we denote spatial averaging in Eulerian 
space, e.g., for a spatial tensor field -A(x, t) — {Aij (x, t)} 
we simply have the Euclidean volume integral normalized 
by the volume of the domain: 



(2) 



For domains T> with constant mass M-p, as for La- 
grangian defined domains, the average density is 



(e(ti)) Vl M z 



(3) 



Averaging Raychaudhuri's equation one finds that the 
evolution of the scale factor a-p obeys the generalized 
expansion lawQ ( Q , see also Appendix |a]) : 



A = Qv, 



(4) 



with Newton's gravitational constant G, the cosmolog- 
ical constant A, and the "backreaction term" Qp. For 
Qt> = this equation equals one of the Fricdmann equa- 
tions for the scale factor a(t) = ap(t) in an homogeneous 
and isotropic universe with uniform density qh = (o)p 
{Qv = is a necessary, but also sufficient condition to 
guarantee ap(t) = a(t)). As we shall see in the next sub- 
section the backreaction term Q-p depends on the inho- 
mogeneities of the velocity field inside V, and is generally 
not zero. 



B. Backreaction in Eulerian coordinates 

Directly following from averaging Raychaudhuri's 
equation we obtain the backreaction Q-p depending on 
the kinematical scalars, the expansion rate 9, the rate of 
shear a, and the rate of vorticity u> (see Appendix [A|): 



Qv 



l{(e 2 )v-(e) 2 v)+H" 2 -° 2 \ 



(5) 



Often it is more convenient to work with the invari- 
ants of the gradient of the velocity field^] as defined in 
Appendix |b|: 



Qv = 2 (IIKj)), 



v ■ 



(6) 



The Eulerian approximation schemes are usually de- 
fined with respect to a reference frame comoving with 
the Hubble flow. To define such a reference frame within 
an inhomogeneous cosmology we have to assume that 



for domains U significantly larger than the averaging do- 
main T> the expansion a(t) is a given solution, e.g., a 
solution of the Friedmann equations. Hence, using co- 
moving coordinates we already assume that Qu = for 
very large domains U. As we will see from Eq. ( |l(i| ) (see 
also Subscct. |VB| ), this can be achieved by setting pe- 
riodic boundary conditions on the inhomogeneous fields 
at the scale of the domain U (see also Q for more pre- 
cision). With the Hubble-parameter H = a/ a we de- 
fine the comoving Eulerian coordinates q := x/a and 
the peculiar- velocity u := v i/x. Using the derivative 
d q -Ui = dui/dqj with respect to comoving coordinates 
we obtain for the first and second invariants: 

IK,-) - 3H + I(ui d ) = 3H + ±1 (d qj u z ) , 
II(v id ) = 3H 2 + 2HI(u id ) + II(u id ) (7) 



3H 2 



The volume-average (A) v of a tensor field A can be 
written as a volume-average over comoving domains T> q 
defined by^ V q — V(t)/a(t) with the comoving volume 
V q = V/a 3 : 



d 3 q A. 



(8) 



The mass is conserved inside this domain V q . A well- 
known example for such comoving mass-conserving do- 
mains are the volumes used in N-body simulations. 

Now, the backreaction term in comoving coordinates 
reads 



Q<D= -TT 



2 -L f d 3 q II (d qjUi ) 



V, 



d 3 q I (9, 



q JT> 



(9) 



i.e., the form of Q-p is such that all Hubble terms in 
the velocity gradient cancel out: only inhomogencitics 
contribute to backreaction. 



Using Eqs. (Bl) and ( JB2| ) we write the backreaction 
as a volume-average over divergences. Hence, using the 
theorem of Gauss we obtain: 



1 



!v = 



2 A- / dS ■ (u(V, • u) - (u • V» 

V Q JdV, 



Vc 



dS-u 



q JdV 



(10) 



1 dt is the total (convective) time derivative, also denoted 
with an over-dot " ' " . 

2 A comma denotes partial derivative with respect to Eule- 3The division V/a is understood point-wise: V/a = 
rian coordinates d/dxi =, i. { x / a I x ^ ^l- 
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with the surface &D q of the comoving domain T> q , the 
surface element dS, and the comoving differential oper- 



ator Vq. From Eq. ( j!0[ ) we directly obtain Q-p = for 
a domain without boundaries, i.e. for periodic peculiar- 
velocity fields (see also the subsection on N-body simu- 
lations |vb| ). 



C. Backreaction in Lagrangian coordinates 

Let f (•, t) : X I— > x denote the mapping which takes the 
initial Lagrangian positions X of fluid elements at time 
ii to their Eulcrian positions x at time t; i.e. x = f (X, t) 
and X := f(X, tj). Then the velocity is v = f(X, t), and 
the Jacobian determinant J of the mapping f connects 
the volume elements d 3 x = Jd 3 X. 

The Lagrangian domain T>\ is connected with the Eu- 
lerian domain by T>\ = f _1 (2?, t), as long as f is unique. 
The spatial average of a tensor field A in Eulcrian and 
Lagrangian space are related in the following way: 



1 1 
Wp~Ji 



d 3 X J A = 



1 



with 



V 
Vi 



( J )c, = 77 = a v- 



{JA) V , (11) 



(12) 



The backreaction can be expressed as a volume-average 
in the initial Lagrangian domain T>\ with the velocity field 
v = f and x = f: 



Q 



V 



(J) 



-{J II^))^--!— (J I(v h] )) Vi 



(13) 



D. A qualitative discussion 



Eqs. (^) and (|5|) show that as soon as inhomogeneities 
are present, they are sources of the equation governing 
the average expansion. Upon integrating Eq. (Q) we ob- 
tain ( Q used a different sign convention for Q-p): 



df, 



Vp 



8irG (g} 1 



'p 



(14) 



kp enters as an integration constant depending on the 
domain T>. With Hp — ap/a-p we define an effec- 
tive Hubble-parameter, and a dimensionless "kinemat- 
ical backreaction parameter" : 



Or 



1 f* 

r— 2- / dt Q-p 2a v a-p 



3 dp-i-ij, j t; 



(15) 



in addition to the common cosmological parameters: 



p 



8ttG 
3l| 



A 



3^' 



1 D 



a p n p 



(16) 



However, here, all -parameters are domain- 
dependent and transformed into fluctuating fields on the 
domain; for Q-p = the cosmic triangle is undistorted 
and the parameters acquire their global standard values. 
Comparing these definitions with Eq. ( |l4] ) we have 







• ng = i. 



(17) 



In Friedmann-Lemaitre cosmologies there is by definition 
no backreaction: fig = 0. In this case a critical universe 
with + = 1 implies fcp = 0, as usual. If instead 
fig ^ 0, and for simplicity f2j? = we have an additional 
kinematical contribution resulting in 1 = + fi^ + f2g. 
Contrary to the standard model, all parameters are im- 
plicit functions of spatial scale: the curvature parameter 
appearing as an integration constant can be different for 
different averaging domains T>, the parameter of the cos- 
mological constant is scale-dependent through H-p. As 
we shall see, f2j? can experience large changes, if initially 
k-p departs (even slightly) from zero. 

Let us consider irrotational flows with uj — 0, which 
we may consider a good approximation until the epoch 
of structure formation. The averaged shear fluctuations 
(ct 2 ).^, > and the fluctuations in the expansion rate 

{6) 2 V -(6 2 ) V = ({9 - (0}p) 2 ) v > compete in the back- 
reaction Eq. (0): 

• Qp is positive if the shear term dominates the ex- 
pansion term. This leads to an accelerated struc- 
ture formation at least in the early evolutionary 
stages. We then may interpret Q-p as a "kinemati- 
cal dark matter" . 

• Qp is negative if the fluctuations in the expansion 
rate dominate. In this case Q-p acts like a (pos- 
itive) "kinematical cosmological constant" leading 
to accelerated expansion. 

Clearly such a "kinematical dark matter" will show a 
time dependence different from ordinary (dark) mat- 
ter; also, if Qp mimicks a cosmological term, its 
time-dependence will also be different than that in a 
A— cosmology. 



III. QUANTIFYING THE BACKREACTION 

In the last section we gave different forms of the back- 
reaction term, and discussed qualitatively the expected 
consequences. Now we quantify them. Exact inhomoge- 
neous solutions for estimating the amount of backreac- 
tion are only available for highly symm etric models like 
for the spherical model (Subsect. Ill A| ), and for a model 
with plane symmetry (Subsect. [II D| ) . In generic situa- 
tions we have to rely on approximations. 
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A. The spherical collapse 

We adopt the usual assumptions. The domain of av- 
eraging T> — Br is a sphere with radius R. The velocity 
v inside Br is only depending on the distance r to the 
origin and always parallel to the radial unit vector e r : 



v(r) e r 



(18) 



Hence, we exclude rotational velocity fields. Moreover 
the domains stay spherical at all times. The averaged 
first invariant may be obtained directly using Gauss' the- 
orem: 



( I ( V i,j))'D = 



AirR 3 



[ dS • v(r)e r = 3 

JdB R 



r ' 



(19) 



whereas the averages of the second and third inv ari ant s 
require some basic calculations. Using Eqs. (B2), (B3) 
and again Gauss' theorem leads to 



1 2 
(ll(v hj )) v =- (I(v id )) v , 

1 3 

(III(v hj )) v =— (I{v hj )) v . 



(20) 



Combining these terms in the backreaction Qb r given by 
Eq. (0) we get for spherical symmetry 



-jsphcrical 



0. 



(21) 



as expected from the Newtonian analogue of Birkhoff 's 
theorem (Newton's "iron sphere" theorem). 



B. Backreaction in the Eulerian linear 
approximation 

We calculate the backreaction in the Eulerian linear 
approximation using Eulerian coordinates comoving with 
the background Hubble flow. In this ansatz we assume 
that there exists a well-defined background density qh 
and a global expansion factor a(t) on very large scales, as 
already mentioned in Subsect. II B. We define the density 
contrast 5 = (g — qh)/ Qh- Additionally we assume that 
the growing mode of the peculiar-velocity field u is dom- 
inating, and we can neglect rotations. The time evolu- 
tion of u(q, t) is then proportional to the initial peculiar- 
velocity U(q) = u(q, t\) (see also Appendix |c]): 



u(q,t) = c(t)U(q) 



(22) 



With the irrotational peculiar-velocity according to 
Eq. (p2|) the invariants may be written as 



l(d qj u l )=c{t)l(d qj U l ) 



II {d, 



\t) II (d q3 Ui 



(23) 



A comoving domain is approximated byQ T> q = 
T>(t)/a(t) w T>i (the Eulerian approximation does not 
strictly conserve mass). In this approximation the vol- 
ume at time t is simply V — a 3 VI, and the backreac- 
tion in the linear approximation can be calculated using 
Eqs. (|), (§, and (f 



Q l V — ~ yQVi 



(24) 



with the backreaction term evaluated on the initial do- 
main 



Q Vi =2 (U(d qj U t )) Vi -- (l(d q] U t )) 2 V; 



(25) 



For an Einstein-de-Sitter background c(t) = (t/ti) 1 / 3 
and a(t) — (t/ii) 2 / 3 , and the backreaction term decays 
in proportion to the global expansion: 



Qv 



-2/3 



Q 



(26) 



Comparing the sources Q£ and {q)^ = QH{U)/ a in 
the generalized Friedmann equation we find growth of 
the deviation from the standard Friedmann models, on 
an Einstein-de-Sitter background: 



Q-D 



With the function 



1 /"* c 2 
h(t) = — dt' — 2a T> a- D , 
o-t, J ti a 



(27) 



(28) 



we can write the backreaction parameter fig, as given in 
Eq. (|l5|), in the following form: 



h(t) 



3 H%, ' 



(29) 



emphasizing the analogy with a time-dependent cosmo- 
logical A-term (see Eq. (|l6|)). 

In addition to using the linear approximation for the 
time evolution of the (irrotational) velocity field we as- 
sume that the comoving domain T> q is only weakly de- 
formed from the initial domain T>\. This approximation, 



4 To make this more rigorous we relate T> q = V(t)/a(t) to an 
initial domain T>i using the "Zel'dovich approximation" f z 
(see Eq. @): 



V q = V(t)/a(t) 
= {q|q = 



i f z (A,t)/a(t) = 
X + e(t)VoV(X)}: 



{q|q = X}=A, 



where the last approximation is valid for I^Vo^l <C 1- Pad- 
manabhan and Subramanian M discuss an extension. 
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applicable in the initial stages of structure formation, is 
not needed in the Lagrangian picture, where we have 
a de finite mapping from T>\ to T>(t) as given in Sub- 
sect. [I C This illustrates that the Eulerian perturbation 
theory is not well-suited to calculate the backreaction of 
inherently mass conserving, i.e. Lagrangian, domains. 



C. Backreaction in the "Zel'dovich approximation" 

Using the Lagrangian perturbation scheme one can 
construct approximate solutions which are able to trace 
the process of structure formation into the Eulerian non- 
linear regime. As a subclass of the first-order Lagrangian 
perturbation series ||] the "Zel'dovich approximation" 
H is frequently used and well-studied both theoretically 
and numerically (see, e.g., |fio|l , and for a tutorial, 
and refs. therein). 

As with the Eulerian perturbation approximation we 
assume a homogeneous-isotropic background model on 
very large scales, specified by the time evolution of the 
global expansion factor a(t). The trajectory field in the 
"Zel'dovich approximation" is given by 



To ease the calculations we define the functional deter- 
minant of three functions ^4(X), -B(X), C(X): 



f z (X,i) 



a(i)(X + £(i)V V(X; 



(30) 



■0(X) is the initial displacement field, Vo the gradient 
with respect to Lagrangian coordinates and £ a global 
time-dependent function (where £ (a) is given for all back- 
ground models a(t) in $2\). 

Given the trajectory field there are two ways of imple- 
menting the approximation for the effective scale factor 
a-p. One is based on the volume deformation of fluid 
elements due to the trajectory field. This results in a 
"passive" estimate of the time-dependence of the volume 
measured by the Jacobian determinant: 

j z (x, t) = a 3 (i + eii + e 2 n ; + e^iii) , (31) 

with the invariants of the initial displacement field^ tyuj'- 
I; := I(V|y), Hi := II(^,), IHi := Hltyfc-)- (32) 
Now (J z ) v . =■ (ap°) 3 is given by 

( a £n)3 = a 3 (1 + £ (I.)^ + ? (11,)^ + e (Hh) Vi ) • 

(33) 

A second possibility is to consider the dynamical back- 
reaction equation itself and approximate the backreac- 
tion term from the velocity field w z = f z using x = i z . 



5 With \i we denote differentiation with respect to the La- 
grangian coordinates X;. 



J{A,B,C) 



d(A,B,C) 
d(X 1 ,X 2 ,X 3 



£i]kA\ t B\jC\ k . 



(34) 



(For example: the Jacobian determinant of the 
Lagrangian mapping f = (/i.,/2,/3) reads J = 
J{fi- /2j /3)-) Now we express the invariants of the ve- 
locity field v = f in terms of functional determinants, 



(35) 



which may be verified using 

V h3 = Jj e jklJ(ft, fk, fi)- (36) 



Inserting the approximation Eq. ( p0[ ) we get with 
b{t) = o(t)£(i): 

J(i Z ,.ff,fk)=aa 2 J{X h X h X k )+bb 2 Ji^i,^,^) 
+ {2aab + a 2 b) J{X i ,X j ^\ k ) 
+ (2abb + ab 2 ) J(X U %,ip ]k ), (37) 

jtffjf, Ik ) = ^ a J{Xi,X j3 X k ) + b 2 b Jty\ u %,i)\ k ) 
+ {2aab + a 2 b) J{X i ,X j ,^\ k ) 
+ {2abb + ab 2 ) JiXi,^,^). (38) 

Again, these functional determinants are related to 
the invariants of the initial displacement field ip^j (see 
Eq. ©): 

^ijkJ{Xi,Xj,X k ) = 6 

e ijk J(X i ,X j ,i) {k )=2I u (39) 

^ijkJ{Xi,i)\j,i)\k) = 2iii, 

£ijkJ(ip\i,ip\j,ip\k) = 6IIIi. 
Hence, with K := £/£ and H := a/a, 

^ eijkJiff, if, /f) =6H + 2 I ; £(3H + K) 

+ 2 II; f(3H + 2K) 

+ 6 III; £ 3 (H + K), (40) 

4 ^kJ{fi , //, fk) =SH 2 + 2 Ii £(3tf 2 + 2HK) 
a 

+ 2 II; £ 2 (3H 2 + AHK + K 2 ) 

+ 6 Illi f(H 2 + 2HK + K 2 ). 

(41) 

Using Eqs. (^, [id]) inserted into Eq. (|35| ) we can calculate 
the invariants of u.y as they appear in Eq. (|l]). Com- 
bined with Eq. ([53]) the backreaction term separates into 
its time-evolution given by £(t) and the spatial depen- 
dence on the initial displacement field given by averages 
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over the invariants Ii , Hi , Illi : 

e 



2 (Hi) 



Pi 



2 (I-Y 

3 \ T>i 



IV-J 



•c (6(111^-1(1^ 



eteih)- (nii) 5 



(42) 



The numerator of the first term is global and corresponds 
to the linear damping factor; in an Einstein-de-Sitter 
universe £ 2 oc a -1 . The denominator of the first term 
is a volume effect, whereas the second term in brackets 
features the initial backreaction as a leading term. 

In the following we will stick to an Einstein-de- 
Sitter background with a{t) = (t/U) 2 / 3 and VoV>(X) = 
3/2U(X)ti |§. £(t) equals a(t) - 1, since f(X,i ; ) = X 
with a(ti) = 1. 

In the early stages of structure formation with <C 1 
we get 



Ql, 



Q-Di, 



(43) 



identical to the early evolution in the Eulerian linear ap- 
proximation Eq. (|2^). Moreover, for t = U we have Q%, — 
Qv ; consistent with Eq. (|J). For late times and posi- 
tive invariants (Ii)^,. , (Hi)x>, , (IHi)x>. > the dominat- 
ing contribution is proportional to £ 2 /£ 4 and Qf, decays. 
However, for negative invariants, the term diverges 
at some time when 1 + £ (Ii)j> + £ 2 (Hi)r>. + £ 3 (IHi)r>. 
passes through zero. This is the signature of pancake- 
collapse on the scale of the domain. A variety of complex 
phenomena may be expected for averaged invariants with 
different signs. 



D. An exact reference solution: the plane collapse 

The "Zel'dovich approximation" is an exact three- 
dimensional solution to the Newtonian dynamics of 
self-gravitating dust-matter for initial conditions with 
II(i/>\ij) = = III(V>|ij) at each point jl]|. This "locally 
one-dimensional" class of motions contains as a sub case 
the globally plane-symmetric solution |Q. Eq. ( |42| ) re- 
duces to 



Qt> 



plane 



2 (Ii) 2 ,, 
'3 (l + £(Ii>p,) 2 



(44) 



We compare the two exact solutions, the plane and 
the spherical model in Subsect. IV A 2. For n egati ve 



, corresponding to over-dense regions (Eq. (Cll)), 

Q!p lanc is diverging at some time when 1 + £(t) (Ii)^,. 
approaches zero. Our special initial conditions imply 
a one-dimensional symmetry of inhomogeneities, and 
the diverging Q!p lanc is supposed to mimick the highly 
anisotropic pancake collapse in the three-dimensional sit- 
uation. 



E. Zel'dovich's approximation and the spherical 
collapse 



In the situation of a spherical collapse (Subsect. [II A) 



we have seen that the averaged invariants obey relations 
resulting in 



(iii)i 



and (Illi) 



1 

27 



(Ii) 



(45) 



Inserting this into the backreaction term Eq. (J42J) we ob- 
tain Q§ R = = <5g P ^ CIlcal . Hence, the "Zel'dovich ap- 
proximation" together with the generalized Friedmann 
equation results also in an exact solution for the evolu- 
tion of the effective scale factor op, if restricted to the 
special initial conditions U = U(r) e r . 



IV. THE EVOLUTION OF COSMOLOGICAL 
PARAMETERS 

The backreaction term itself decays in the linear ap- 
proximation and also at early stages in the "Zel'dovich 
approximation" . However, to quantify the deviations 
from the behavior of the scale factor in the standard 
model, the different strength of the sources in the gen- 
eralized Friedmann equation has to be compared. As 
already noted in the Eulerian linear approximation, the 
dimensionless contribution to backreaction (as compared 
to the global mean density) grows. 



A. Evolution of the scale factor ax> 

We proceed as follows. We use the backreaction terms 
calculated in the preceding sections, insert them into the 
generalized expansion law Eq. (Q), and solve for a-p(t). 
We determine generic initial conditions numerically as 
summarized in Appendix ^|. 

Assuming an Einstein-de-Sitter background we sub- 
tract a standard Friedmann equation 3^ + A-kGqh = 0, 

from Eq. (^), use qh = §^77 and obtain the following 



differential equation for a-pit): 



3( av_a 



(46) 



with (5} v = {{q)x, — Qh)/qh specifying the averaged 
density contrast in T>. For Q-p = and (8)^ = this 
equation simply states that the time evolution of a do- 
main follows the global expansion, ax>{t) — a(t). The 
averaged density contrast is relat ed t o the averaged first 
invariant as given in Eq. ( |5ll| ). For Q v = and 

1 + (fy-n = the evolution of a-p is still of Fried- 

mann type, but with a mass different from the homo- 
geneous background mass. An important sub case with 
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Qt> = and (6) v ^ is the well-known spherical top-hat 
model Jl5|. In Eq. ( ^6|) there are two sources determin- 
ing the deviations from the Friedmann acceleration, the 
over/under-density and the backreaction term: the evo- 
lution of a spatial domain is triggered by an over/under- 
density and velocity fluctuations. Eq. ( ]46| ) can be inter- 
preted as a natural generalization of the top-hat model. 
Engineer et al. |l6| also incorporate shear and angular 
momentum to formulate and study generalizations of the 
spherical top-hat model. Their emphasis, however, is fo- 
cussed on the improvement of the virialization and other 
clustering arguments. 

To solve this differential equation we specify the initial 
conditions according to the background model: 



a(t) 



2/3 



with ti 



3H a (l + Zi 



|3/2 - 



(47) 



Hq is the present value of the Hubble-parameter of the 
background model, and z\ the redshift where we start 
our calculations. With Ho = 50km s _1 Mpc~ 1 and z\ = 
200 the initial starting time is U = 4.6 10 _3 Gy. The 
scale factor a%> of the domain is chosen to match the 
background scale factor a at time ti: 



av(ti) = a(ti) = 1. 



(48) 



We stop at the present epoch z = and t n = 13Gy, 
with a(to) = 201. Following directly from Eq. (A3), a-p 
is given by 

av{U) = d(t,) (1 + § = |: (1 + 3 (li)^) • (49) 

Often ax>{ti) = a(ti) is used to specify the initial con- 
ditions for the investigation of collapse times within the 
spherical model. Numerical tests showed that for these 
inconsistent initial conditions the collapse is delayed as 
already discussed by Bartelmann et al. fl7|| . 

For the numerical integration we used the Runge- 
Kutta routine given in Press et al. |Q . Identical results 
were obtained with the NDSolve routine of Mathemat- 

ICA. 



The ensemble mean E, i.e. the expectation value, of the 
volume-averaged invariants vanishes in conformity with 
our global set-up: 

E [(IObJ = E [(IIObJ = E [<HIi) e J = 0. (50) 

However, for a specific domain, any of the volume- 
averaged invariants may be positive or negative. These 
invariants fluctuate with the variance, e.g., Cj(i?) = 



e.g. 

. We consider one-a fluctuations, e.g. (Ii) 



B, 



±di(R), and also two- and three-tr fluctuations of the 
averaged invariants in our calculations of the time evolu- 
tion of ax>(t). Clearly, these fluctuations depend on the 
shape of the power spectrum and the radius of the ini- 
tial domain. In Table. | we give the values for spherical 
domains of different sizes, as calculated in Appendix ^ 
for the standard CDM model. Moreover, we could show 
that the volume-averaged invariants are uncorrelated for 
a Gaussian initial density field, and may be specified in- 
dependently. 

TABLE I. The mean mass fluctuations ai(R), au(R), and 
oin(-R) for an initial domain of radius R are given, calculated 
using the CDM power spectrum. To make these numbers 
more accessible, also the linearly extrapolated mean fluctu- 
ations a(tri)ai(R) for a domain with scaled radius a(to)R at 
present time are given (compare Fig. [I]). 

aR [Mpc] 



16 50 100 251 



503 



1005 



aai 



2.6 1.0 0.27 0.1 0.025 0.0068 0.0019 



R [Mpc] 0.025 0.08 0.25 0.5 

(Ji xlO 3 13 5.0 

cm xlO 6 90 20 

am xlO 12 81000 4600 



1.4 0.51 
2 0.4 
100 5 



1.25 
0.13 
0.07 
0.08 



2.5 5 

0.03 0.01 

0.02 0.004 

0.001 0.00004 



1. Initial conditions 



In Sect. Ill we expressed the term Qx> using globally 
given time-functions a(t), £(i), and using the volume- 
averaged invariants (Ii)-p. , (Ili)^,. , (Illi)^. for the initial 
domain D\. The relation of these volume-averaged in- 
variants to the initial power spectrum is discussed in Ap- 
pendix |C| in detail using the linear approximation. We 
determine the initial values of the averaged invariants 
for a spherical initial domain Br of radius i?, assuming a 
Gaussian initial density field, with a Cold-Dark-Matter 
(CDM) power spectrum. 







2. The plane and spherical collapse 

We calculate the evolution of the scale factor a-p(t) 
according to Eq. (^) using the backreaction term cal- 
culated within the "Zel'dovich approximation" ( fi~2|) with 
a(t) = (t/U) 2 / 3 and £(t) = a(t) - 1. 

First we consider initial conditions with = = 

and 7^ 0, the plane collapse solution 

Eq. (fl4|). We choose comparable to the mean 

fluctuations o\ (R) in initially spherical domains with ra- 
dius R £ {0.025, 0.08, 0.25, 0.5}Mpc (a scaled, present- 
day radius of a(t )R 6 {5, 16, 50, 100}Mpc) as given in 
Table H. From Fig. || it can be seen that the evolution 
of the scale factor a-p {t) may significantly differ from the 
background expansion a(t) for an initial domain radius 
R = 0.25Mpc (a scaled radius of a(t )R = 50Mpc). For 
small domains with larger (Ii) B the effect is more pro- 
nounced. For negative (an over-dense region) the 
region shows a pancake collapse with Q-p — > oo, as a re- 
sult of the one-dimensional symmetry in the initial condi- 
tions. For initial domains with R > 2.5Mpc (a scaled ra- 
dius of a(t )R > 503Mpc) the value of (Ii)e R is typically 
so small that the ax> {t) closely follows the Friedmann so- 
lution until present. In the following we compare the two 
inhomogeneous exact solutions, the plane and spherical 
collapse, disentangling the influence of the backreaction 
term Qx> from the influence of the over-density. 



2 r -;0 



200 - 



150 - 




100 - 



t [Gy] 

FIG. 2. The curves show the scale factor ax>(t) in the 
spherical collapse (dotted) and the plane collapse (solid) for 
(Ii) BR e {-13,-5.0,-1.4,-0.51,0.51,1.4,5.0,13} x 10~ 3 , 
bending up successively. These values correspond to one-cr 
fluctuations of spherical domains with a scaled radius of 
a(t )R £ {5, 16, 50, 100}Mpc (see Table |). The homogeneous 
background solution is given as the dashed line. 



zero for a spherical domain with spherical symmetric ve- 
locity field. However, this domain may still be over- or 
under-dense compared to the background density. The 
evolution of the scale factor of this domain is then deter- 



mined by Eq. (|q) with Qx> = and 



which is 



exactly the spherical collapse model [|15| 

In Fig. we compare the evolution of a,T>(t) in this 
spherical collapse model with the plane collapse model. 
For | I = 1-4, typical for domains with a scaled ra- 

dius of a(to)R = 50Mpc, both models give nearly the 
same evolution. However for increasingly over-dense re- 
gions (negative (Ii) 0R ) the plane model decouples and 
collapses significantly earlier than the spherical model. 
Also the expansion of under-dense regions is slowing 
down in the plane compared with the spherical collapse 
model. These results are in agreement with earlier re- 
sults on collapse-times of structure in Lagrangian per- 
turbation schemes 



It 



22 1 



3. Generic initial data 

The behavior of ax> (i) becomes more interesting if we 
consider the other two invariants, going beyond the ex- 
act solutions. For simplicity we choose (Ii) B = corre- 
sponding to a domain with {S) B = matching the back- 
ground density at initial time t{. The effect of a nonzero 
second invariant (IIi) B ^ is visible in Fig. ||. Al- 
ready in a domain with an initial radius of R = 0.25Mpc 
(a scaled radius of a(to)R ~ 50Mpc) the a-p may devi- 
ate from the Friedmann solution. Although we consider 
volumes with no initial over-density and with only one- 
cr fluctuations of (Hi)g , domains of a current size of 
a(to)R~ 16Mpc may collapse. 
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As shown in Subsect. Ill A, the backreaction Qt> is 
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FIG. 3. The curves shown are the evolution of the scale 
factor axj(i) for initial conditions with (Ii) B/? = = (IIIi) Bfl 
and (Il^g G {-90,-20,-2,2,20,90} x 10~ 6 shown as the 
solid lines bending up successively. These values correspond 
to one-cr fluctuations of spherical domains with a scaled ra- 
dius of {5,16,50}Mpc (see Table |l|). The homogeneous back- 
ground solution is given as the dashed line. 

The influence of the third invariant on the evolution 
of the scale factor is small (in the temporal range con- 
sidered) compared to the influence of (Ii),g and (IIj) g . 
To illustrate this we consider domains with (Ii)g = = 
(Ilj)^ and one-cr fluctuations of (IIIj) Br . For an initial 
domain of radius 0.025Mpc (a scaled radius of 5Mpc) 
a significant influence is visible (Fig. ||). Already for 
an initial domain of radius 0.08Mpc (a scaled radius of 
16Mpc) the influence is small, becoming negligible on 
larger scales. 
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FIG. 4. The evolution of the scale factor ap(t) 

for initial conditions with (Il)g R = = (IIi) gR and 
(IIIi) B G {-81,-4.6,4.6,81} x 10" 15 shown as the solid 
lines bending up successively. These values correspond to 
one-cr fluctuations of spherical domains with a present day 
scaled radius of 5 and 16Mpc (see Table Q). The homoge- 
neous background solution is given as the dashed line. 

Up to now we only considered domains with initial 
conditions chosen as one-cr fluctuations in either (Ii)^ , 
(Hi) BR , or (IIIj)g R . As can be seen from Fig. || the 
evolution of the scale factor is strongly depending on the 
initial conditions, even for a domain with initial radius 
0.25Mpc (a scaled radius of 50Mpc), if we also consider 
two-cr or three-cr fluctuations. Such a domain may even 
collapse. 
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FIG. 5. The evolution of the scale factor oz>(t) for spher- 
ical domains with a scaled radius of 50Mpc. We consid- 



ered one-cr (dotted line) 
three-cr (long dashed line 
and (IHi) B (see Table 
tained for (h) BR , (IIi) B 

(Ii> s ,(IIi) B ,(IIIi> R „ ; 



two-cr (short dashed line), and 
fluctuations of (1x} Br , (IIi) b , 
. The evolution shown are ob- 
(IIIi) e < (over-dense) and 
B R > \"i/8 R ' V^'/Bfl ^ (under-dense). The homoge- 
neous background solution is given by the solid line. 



Using the present approach we already assume that 
the backreaction term gives a negligible contribution on 
the largest scale. Hence, for increasing size of the initial 
domain the evolution of the scale factor goes more and 
more conform with the Friedmann evolution (see Fig. |^). 
Still for a domain with a scaled radius of lOOMpc devi- 
ations of the order of 15% in the scale factor at present 
time may be obtained for initial fluctuations at the three- 
cr level. The scaled volume of such a domain corresponds 
to a cube with w 160Mpc side length - still a typical 
volume used in A-body simulations. For domains with 
a scaled radius of 250Mpc (« 400Mpc side length of a 
cube) we may get deviations of the order of 3% from 
the Friedmann value a fa), for init ial flu ctuations at the 
three-cr level. In Sects. IV A 5| and IV B we will see that 
the influence on the domains' Hubble- and deceleration- 
parameter as well as on its density parameters is more 
pronounced, because these parameters invoke the time- 
derivatives of ax>(t). 
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FIG. 6. The evolution of the scale factor av>{t) for spheri- 
cal domains with a scaled radius of lOOMpc (upper plot) and 
250Mpc (lower plot). The same conventions as in Fig. H are 
used. 



the "Zel'dovich approximation" is reasonable for domains 
with a scaled radius larger than lOOMpc. The Eulerian 
linear approximation fails to describe a collapsing do- 
main. However, we see that the backreaction term is not 
only important for small domains, where the nonlinear 
evolution plays a dominant role, but also for large do- 
mains within the realm of linear theory. We emphasize 
that these remarks are concerned with the average per- 
formance; linear theory does not predict correctly where 
mass is moved, as cross-correlation tests with Lagrangian 
perturbativc schemes and N-body runs (even if smoothed 
on larger scales) have shown |23| . 




5 10 
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FIG. 7. The evolution of the scale factor ax>(t) for three-cr 
initial conditions as already used in Fig. ^j. The curves are for 
over-dense domains with a scaled radius of 50, 100, 250Mpc, 
and under-dense domains with a scaled radius of 250, 100, 
50Mpc, bending up successively. The evolution of af> is given 
by the solid lines, the evolution of a^, n by the dotted lines. 
The evolution according to the Friedmann solution is given 
by the dashed line. 



5. Evolution of the Hubble- and deceleration-parameter 



4- Comparison with the Eulerian linear approximation 

Using the backreaction term ( p6| ) in the generalized 
evolution equation (jig) we calculate the time evolu- 
tion of ai^(t) in the Eulerian linear approximation. In 
Fig. ^ we compare the evolution of the scale factor in 
the "Zel'dovich approximation" with the evolution in the 
Eulerian linear approximation. Even though we consider 
three-cr initial conditions, the evolution of under-dense 
domains with a scaled radius larger than 50Mpc is de- 
scribed nearly perfectly by the Eulerian linear approx- 
imation. For over-dense domains the concordance with 



Not only the evolution of the scale factor itself is inter- 
esting, but also the derived quantities like the Hubble- 
parameter H-p = d-p/a-D and the deceleration parameter 
qv = ~^7tt^~ are of specific interest. 

r-i 

From Fig. |8| we see that at present a 20% variability of 
the Hubble-parameter ifp is quite common in a domain 
with a scaled radius of lOOMpc (three-cr initial condi- 
tions). For a domain with a scaled radius of 250Mpc still 
3% variation is possible. This agrees well with the results 
obtained by Wu et al. 1 24| for a regional low-density uni- 
verse (see also |E5|, B6], and refs. therein). 
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FIG. 9. The evolution of the deceleration parameter 
qv(t) for a spherical domain with a scaled radius of lOOMpc. 
Over-dense domains show an increased qv compared to an 
Einstein-de-Sitter universe (solid line) , under-dense domains 
a decreased qv- The same conventions as in Fig. |H] apply. 

Consider the situation that our regional Universe may 
be an under-dense portion (say, resulting from an under- 
density at 0.5-cr on the scale of the considered portion). 
Then, we may well live in a shear-dominated part of 
the Universe: many low-amplitude anisotropic struc- 
tures surround us which are compatible with a regional 
under-density, but would suggest a dominance by shear 
fluctuations (say, resulting from a positive second invari- 
ant in the initial conditions at three-cr). Such domains 
with density-contrast and velocity fluctuations working 
"against" each other may show a qualitative change in 
their evolution as illustrated in Fig. [llj: an initially 
under-dense region, however, dominated by shear fluc- 
tuations may show a deceleration. 



FIG. 8. The evolution of the Hubble-parameter Hx>{t) 
for a spherical domain with a scaled radius of lOOMpc. 
Over-dense domains show a smaller Ht> than in an Ein- 
stein-de-Sitter universe (solid line), under-dense an increased 
Hx>- The same conventions as in Fig. |^ are used. 

The effect becomes more pronounced if one considers 
the deceleration parameter qp, as shown in Fig. [)| In 
a domain with a scaled radius of lOOMpc the q-p may 
show a present day value of more than 200% different 
from the background Einstein-de-Sitter universe. Even 
for a domain with a scaled radius of 250Mpc a difference 
larger than 15% is possible. 
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FIG. 10. The evolution of the regional deceleration pa- 
rameter <Jx>(t) for a spherical domain with a scaled radius of 
lOOMpc. The initial conditions are chosen as 0.5-cr for 
and three-cr for both and (IIIi) g with the opposite 

sign of (Ii) B . The initially accelerating, i.e. under-dense re- 
gion > 0) finally decelerates (dotted line) whereas the 
initially over-dense region (dashed line) shows an accelerated 
expansion at late times. 



B. Evolution of the density parameters 

Now we turn our attention to the domain-dependent 
"cosmological parameters" 



defined in Subsect. [ID. 
Knowing the time evolution of a-p for given initial condi- 
tions we can calculate the time evolution of 0^, fij^, and 
ilq. Again we assume that the global evolution follows 
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an Einstein-de-Sitter model with fl\ = = f2? and on 
the largest scale flk = 0. However, as shown below, for 
finite domains T> the time evolution of the "cosmological 
parameters" can substantially differ from the Friedmann 
evolution even for small, seemingly negligible backreac- 
tion term. We have for the density and "curvature" pa- 
rameter: 



av(t) 2 a v {t) 

k-D 



av{t) 2 



(51) 



l)ap(ti) 



The evolution of fig is 



with k v = (fi£(i £ ) 
given by Eq. (|i~5|). 

For a homogeneous and isotropic distribution of mat- 
ter Q]? may be related to the total energy inside the do- 
main T> within the Newtonian interpretation. This is not 
possible in the more general inhomogeneous setting we 
considered. Now k-p is an integration constant given by 
the initial conditions. The role of fcp and its relation to 
the average curvature in the relativistic framework will 
be explained in Sect. VA. 



A remarkable result is that for an accelerating i.e. 
under-dense region, f2g may be negligible as seen in 
Fig. [ll], but dramatic changes in the other parameters 
are observed. For a collapsing domain with a scaled ra- 
dius of lOOMpc the mass parameter of the domain may 
even differ by more than 100% from the global mass pa- 
rameter (Fig. 02). 
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FIG. 11. The evolution of the "cosmological parameters" 
(solid line), Of? (short dashed), and Q,q (long dashed) 
in an expanding (under-dense) domain with initial radius 
0.5Mpc (a scaled radius of lOOMpc). The dotted line is mark- 
ing Qj^+flJ? + Qq- The left plot is for one-cr, the right plot is 
for three-o" fluctuations. The corresponding evolution of the 
scale factors is shown in Fig. H. 
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FIG. 12. The same quantities as in Fig. |Tl| are shown, 
now for a collapsing domain. 



C. Self— consistency of the approximations 

In the preceding subsections we calculated the evolu- 
tion of the scale factor a-u by integrating the generalized 
Friedmann equation (|J) for a given backreaction term 
Qv- We used exact backreaction terms for the spheri- 
cal collapse and the plane collapse, but also the approxi- 
mate backreaction term Q|j for generic initial conditions. 
This last approach must not lead to consistent results. In 
the following we test for self-consistency in the spirit of 
Doroshkevich et al. J27j and show that our results ob- 
tained with the "Zel'dovich approximation" are reliable. 

Using and integrating Eq. ( [i"6| ) we obtain an ap- 
proximate solution a% for the scale factor. There is 
also the "passive" volume based approximation for the 
scale factor in Zel'dovich's approximation according to 
Eq. (H|). This kinematical a!g n does not describe the 
same evolution Both solutions are exact for the 

plane collapse initial conditions Hi = = III,. More- 
over, a| reproduces the (exact) spherical collapse (see 

kin 
V 



Subsect. Ill E| ), whereas the kinematical solution a 
does not. Hence, we typically rely on 
The normalized difference 



4 a \ 



(52) 



may serve as an internal consistency check. For initial 
domains larger than 0.08Mpc (a domain with scaled ra- 
dius of 16Mpc) the ei stayed well below 0.1 at all times, 
as long did not approach zero. 

If a| differs significantly from the Friedmann solution 
a we are also interested in the impact of the inconsistency 
between the two solutions a^, and a!g n on the actual ob- 
served deviation: 



-,kin I 



C-2 



(53) 



For all situations with |af, — a\ 3> 0, and especially for 
the cases considered in Fig. £2 is always smaller than 
0.15. 

Hence, the results presented above are self-consistently 
correct within the Zel'dovich approximation. Moreover, 
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we can reproduce exact solutions for both the spherical 
and the plane collapse, i.e. for two orthogonal symmetry 
requirements. 



V. SOME REMARKS 

A. Clues from general relativity 

The advent of the "backreaction problem" has its roots 
in general relativity [Q. Certainly a general relativis- 
tic treatment has to be considered when talking about 
the global effect of backreaction. We have shown that a 
Newtonian treatment falls short of explaining the global 
effect. We remark that the same applies for a general rel- 
ativistic treatment that is designed as close to the New- 
tonian one as possible. For example, Russ et al. ]2g| ] 
have worked with perturbation schemes on a 3-Ricci flat 
background employing periodic boundary conditions on 
some large scale in order to estimate the global backreac- 
tion effect. In their investigation the backreaction term 
(Eq. ( j57l) below) vanishes by construction on the period- 
icity scale and its global effect cannot be estimated (see 
H for further discussion). The question whether an av- 
eraged space section can be identified at all with a Fried- 
mann-Lemaitre cosmology on some large scale, even ap- 
proximately is open and can probably be answered only 
on the basis of a full background-free relativistic investi- 
gation. 

Although it is of particular importance to find strate- 
gies for the estimation of this global effect, the impact 
of backreaction on small scales is also interesting from a 
general relativistic point of view. We will discuss hints 
in this subsection that show that, even on scales that 
are previously thought to be satisfactorily described by 
Newton's theory, general relativity can, and we think, 
will play an exciting role. On top of previous estimations 
of the global backreaction effect in GR we have to ac- 
knowledge that the aspect added by Buchert and Ehlcrs 
Q , namely that backreaction terms appear naturally also 
by averaging out inhomogeneous matter distributions in 
Newtonian theory, is paraphrased in the general relativis- 
tic framework as well. We shall concentrate below on this 
correspondence. 

Let us first briefly review the corresponding equations 
that rule the averaged dynamics in general relativity || . 
Given a foliation of spacctime into flow-orthogonal hy- 
persurfaces we can derive equations analogous to the gen- 
eralized Friedmann equations discussed in this paper (re- 
stricting attention to irrotational flows). Spatial averag- 
ing of any scalar field "J is a covariant operation given a 
foliation of spacetime and is defined as follows: 



— / Jd 3 X ^{t,X l 
Vv Jv 



(54) 



with J := v/ det(<7ij), where gij is the metric of the spatial 
hypersurfaces, and X 1 are coordinates that are constant 



along flow lines. The spatially averaged equations for the 
scale factor op, respecting mass conservation, read: 
averaged Raychaudhuri equation: 



Via 3 



A = Qv; 



v 



averaged Hamiltonian constraint: 



av 



8ttG M v (K) q 



3 Vi<4, 6 



A 
~3 



Qv 
6 • 



(55) 



(56) 



where the mass Mv, the averaged spatial Ricci scalar 
(R-)v and the "backreaction term" Qv are domain- 
dependent and, except the mass, time-dependent func- 
tions. The backreaction source term is given by 



Qv := 2 (!!)„ - | (lf v = | 



(57) 



Here, I and II denote the invariants of the extrinsic cur- 
vature tensor that correspond to the kincmatical invari- 
ants we employed. The same expression (except for the 
vorticity) as in Eq. (||) follows by introducing the split 
of the extrinsic curvature into the kinematical variables 
shear and expansion (second equality above). 

We appreciate an intimate correspondence of the GR 
equations with their Newtonian counterparts (Eq. (^) 
and Eq. (14)). The first equation is formally identical 
to the Newtonian one, while the second delivers an ad- 
ditional relation between the averaged curvature and the 
backreaction term that has no Newtonian analogue. This 
implies an important difference that becomes manifest 
by looking at the time-derivative of Eq. ( |56| ) . The inte- 
grability condition that this time-derivative agrees with 
Eq. (pq) is nontrivial in the GR context and reads: 

dtQ-D + 6— Qv + d t (TZ) V + 2— (K) v = 0. (58) 
av av 

The correspondence between the Newtonian fcu-para- 
meter with the averaged 3-Ricci curvature is more in- 
volved in the presence of a backreaction term: 

(59) 

The time-derivative of Eq. ( |5^ ) is equivalent to the 
integrability condition Eq. ([58]). Eq. ([58]) shows that av- 
eraged curvature and backreaction term are directly cou- 
pled unlike in the Newtonian case, where the domain- 
dependent fc-p-parameter is fixed by the initial condi- 
tions. For initially vanishing kv in Newton's theory, the 
"curvature parameter" will always stay zero. This is not 
the case in the GR context, where backreaction produces 
averaged curvature in the coarse of structure formation, 
even in the case of domains that are on average flat ini- 
tially. In view of our results that the fcp-parameter could 
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play an important role to compensate for under-densities 
in some domain in a globally flat universe, the GR con- 
text suggests that we shall have a more complex situation 
with regard to the backreaction term itself. It remains 
to be seen whether this coupling increases the quantita- 
tive relevance of the backreaction term itself, and how it 
affects the dynamics of spatial domains also on "Newto- 
nian scales". This will be the subject of a forthcoming 
work. 



B. N— body simulations 

Most cosmological N-body simulations solve the New- 
tonian equations of motion in comoving coordinates typ- 
ically inside a cuboid region C with periodic boundary 
conditions. Hence, the boundary of C is empty, dC = 0, 
and from Eq. (|l(]) we directly obtain Qc — 0. For the 
definition of comoving coordinates a background expan- 
sion factor a(t) is required. Consistent with Qc = the 
time evolution of ac (t) is assumed to follow the time evo- 
lution of a Friedmann-Lemaitre background model a(t) 
throughout the whole simulation process. Hence, N-body 
simulations give the correct results (within their resolu- 
tion limits) for the evolution of inhomogeneities in do- 
mains where no backreaction is present at all times. As 
we have seen in the preceding sections, for generic ini- 
tial conditions the backreaction may influence the global 
expansion of domains even with a present size of hun- 
dreds of Mpc's. The difference in the time-evolution of 
a-p compared to a acts as a time dependent source for 
the peculiar-potential </>: 

A„0 = AttGo? qh5 + 3 — (a-pa — a-pa) . (60) 

A general domain will also change its shape as illustrated 
in Fig. |l|. This will also change the evolution of inhomo- 
geneities inside these volumes. Therefore, N-body simu- 
lations only trace a limited subset of the full space of solu- 
tions. Moreover, the fluctuations in the initial conditions 
are und erest imated for periodic boundaries as illustrated 
in Sect. C 3. 



It will be a challenge to incorporate backreaction ef- 
fects in N-body simulations. In order to improve the 
estimation of the backreaction effect, one could embed 
high-resolution N-body simulations into a model for the 
large scales such as the model investigated here. In this 
line Takada and Futamase [^9| have suggested a hybrid 
model that incorporates the full nonlincarities on small 
scales, while the large scales are described by perturba- 
tion theory. A similar procedure was already used for 
high-resolution cluster simulations, surrounded by in- 
homogeneities determined from a low-resolution cosmo- 
logical simulation [^0|. However, abandoning the Fast- 
Fourier-Transform with its periodic boundaries, as typ- 
ically used for solving the Poisson equation in PM- and 
P 3 M-simulations, will be only one technical obstacle. In 



this respect tree-codes are more flexible (see e.g. ||T|] for 
the treatment of periodic boundaries). 

Another question is, whether we are able to estimate 
the effect of the backreaction using simulations. Con- 
sider a partitioning T> n C C of the periodic box C into 
N subsets intersecting only at their boundaries with 
C = \Jn—iT> n . Typically these V n are sub-cubes of the 
periodic box. We obtain for the full backreaction from 
Eqs. © or (§): 



A' 



\C\Q C = Yl \V n \Qv n = 0. 



(61) 



n=l 



In general each Qv n 0, but they cancel over the whole 
simulation box. Let us assume that all the V n have the 
same shape and volume \T> n \ = |C|/AT. According to 
Eq. (^l|) the mean backreaction for such domains esti- 
mated by 



1 N 



(62) 



will always be zero in simulations. An estimate of the 
variance is possible by 



1 N 

71=1 



(63) 



However, the backreaction terms Qv n of the domains 
V n are correlated through Eq. (|6l|), and Eq. ( |63| ) under- 
estimates the variance of the backreaction especially for 
domains T> n comparable in size to C. 



VI. DISCUSSION, CONCLUSIONS AND 
OUTLOOK 

We presented a quantitative investigation of the back- 
reaction effect, i.e., the deviations of dynamical proper- 
ties of an averaged inhomogeneous model from that of a 
homogeneous-isotropic standard cosmology. We did this 
for scalar observables on finite spatial domains such as 
the scale factor and its derivatives (the expansion and de- 
celeration rates), as well as for derived parameters such 
as the cosmological density parameters defined on the 
domains of averaging. 

Although we did this for generic initial conditions com- 
mon in the modeling of large-scale structure, our ap- 
proach is limited by the Newtonian approximation, the 
setting of periodic boundary conditions on the inhomo- 
geneities on some large scale, and the restriction to a 
first-order Lagrangian perturbation scheme for the mod- 
eling of inhomogeneities. We may consider the proposed 
approximation as a prototype model of any averaged in- 
homogeneous cosmology asking for refinements in the re- 
strictions mentioned. 
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First of all we learned that the dynamics of a domain 
is controlled by two effective sources for the expansion 
rate governed by a generalized Friedmann equation: the 
first is the over-/under-density compared with the global 
density, usually studied by the spherically symmetric 
top-hat model. This effect can be absorbed into the 
initial conditions and the average model remains within 
the class of the standard Friedmann-Lemaitre solutions. 
The second is the backreaction term itself that measures 
the departure from the class of standard models. It con- 
sists of fluctuations in the kinematical parts of the veloc- 
ity gradient, the expansion rate and the rates of shear 
and vorticity. Consistently, the spherical model with 
an over- or under-density inside the spherical domain 
describes the isotropic deviations from the Friedmann- 
Lemaitre evolution, as expressed by the vanishing back- 
reaction term. This interpretation is valid if one con- 
siders the evolution of the whole spherical domain. Lo- 
cally however, we have contributions from the expansion 
rate 9 and the shear cr, canceling only on average in the 
spherical domain. Additionally to these isotropic devia- 
tions the plane-collapse or the more general models allow 
for anisotropic deviations from the Friedmann-Lemaitre 
evolution. One may express this directly in terms of the 
averaged invariants: (I)^ = 3H-p is the isotropic devia- 
tion from the Hubble expansion 3H; (II) p incorporates 
the shear responsible for anisotropic deviations. Both 
the fluctuations in the expansion rate (I)p and (II) v 
contribute to the backreaction term. 

On the basis of a first-order Lagrangian perturbation 
approach we proposed an average model for generic ini- 
tial conditions that is exact in two orthogonal cases, the 
plane and the spherically symmetric evolution. Accord- 
ing to these properties the new model generalizes the 
spherical top-hat model. Comparisons showed that the 
collapse is typically accelerated in this more general set- 
ting. At least down to the limit scale of the evolution 
model our results should be correct, i.e. for truncated 
high-frequency end in the initial power spectrum this is 
roughly the scale of galaxy clusters p^] . 

By evaluating the expectation values for the fluctua- 
tions we quantified the expected magnitudes of the parts 
of the backreaction term given a power spectrum of initial 
density fluctuations (in our case we have chosen a stan- 
dard CDM power spectrum). The effect may be viewed 
from two perspectives: a one-cr fluctuation of a given part 
of the backreaction on some smaller scale may be alterna- 
tively viewed as a two-cr fluctuation of a correspondingly 
larger domain. E.g., a negative one-cr fluctuation in all 
parts on a domain of a current radius of lOOMpc (with 
a normalized Hubble-parameter of 0.5, a global density 
parameter of unity and a zero global fc— parameter) will 
give today a matter density parameter of 0.84 compen- 
sated by a k— parameter of roughly 0.15. We see in this 
example that the backreaction parameter itself can be 
negligible quantitatively (fig = 0.01), but qualitatively 
it has a strong impact on the change of the other pa- 
rameters. We emphasize that this is a large effect, since 



a 15% change of the standard parameters on a corre- 
sponding cubic volume of about 160Mpc is a lot for the 
"innocent" expectation values chosen at one-cr. 

These initial fluctuations driving the non-standard 
evolution may arise from modes of the density fluctua- 
tions inside or outside the domain, their relative strength 
depending on the size of the domain and the shape of 
the power-spectrum considered. Using the CDM power- 
spectrum on small scales, where modes larger than the 
domain with a radius of lOOMpc are set to zero, the fluc- 
tuations of e.g. the first invariant are reduced by a factor 
of two to three (see Fig. [l3]). Clearly the large-scale 
fluctuations give the major contribution, but also fluctu- 
ations inside contribute at least 30% in this case. Hence, 
the feasibility of a split between small- and large-scale 
contributions critically depends on the spectrum and the 
size of the volume considered. 

Consider, e.g., the Marklll and SFI catalogs of 
peculiar-velocities. Adopting our Hubble constant these 
samples roughly correspond to a radial depth only 
slightly larger than lOOMpc. According to what has 
been said before we may likely live in an "untypical" 
region when talking about lOOMpc. (The likelihood of 
such events has not been quantified in the present pa- 
per, but is a subject of ongoing work). We therefore 
think that matching results from small-scale data with, 
e.g., the supernovae constraints (see e.g. |53|) has to be 
premature unless we do not have either a considerably 
larger sample, or else, know the fluctuation properties of 
the regional Universe by other observations. Since the 
value of the backreaction term Qx> depends on the veloc- 
ity field inside £>, these peculiar-velocity catalogs may 
offer the possibility of estimating Qx>. In this line two 
papers are of particular interest: by taking the sampling 
anisotropics of the velocity field explicitly into account, 
Regos and Szalay |Q found a large effect (40%) of the 
dipol and quadrupol anisotropics on the estimated bulk 



351 ; using the Eu- 
61 showed that the 



flow of the elliptical galaxy sample o_f 
lerian linear approximation Gorski 
velocity field is significantly correlated even on scales of 
lOOMpc. 

Furthermore, we may well deal with an underestimate 
of the effect, since our assumptions are conservative in a 
number of respects: 

The "Zel'dovich approximation" reliably traces large- 
scale features but lacks structure on small scales which 
would add up in the averages, since shear fluctuations 
and expansion fluctuations are positive definite. Simi- 
lar remarks apply to the small-scale end of the power 
spectrum, which differs between models. Certainly, the 
spectral index on large scales and the position of the 
bend-over plays a crucial role. Thus, as a first point 
we may say that in the family of models usually studied 
the standard CDM model evolved with the "Zel'dovich 
approximation" may be regarded as conservative in these 
respects. 

We started our calculations at z\ = 200 since at higher 
redshifts the (linear) power-spectrum is still changing its 
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shape. Therefore, we have to assume that the evolution 
of the domain follows the Friedmann equation with no 
backreaction up to z\ = 200. Our matter model (dust) is 
not a good description for earlier times, however extrapo- 
lating and starting at z\ = 1000, the resulting deviations 
increase only slightly, well below the percent range for 
the examples considered. 

As a third point, a conservative constraint that we have 
put in is the existence of a Hubble flow on some large 
scale. From our discussion of the general relativistic case 
it is suggested that a realistic inhomogeneous cosmology 
will take its freedom to drift away from the standard 
model. Assuming Qx> — for all times means that we 
force the effect to add up to zero within the whole model 
for all times. This certainly hints towards an increase of 
the effect under study in more realistic situations. 
Fourth, due to the Newtonian investigation we were also 
conservative: in general relativity, the curvature param- 
eter is directly coupled to the backreaction parameter. 
Since we have learned that the Newtonian fc-p-parameter 
plays a crucial role, this remark means that we may 
underestimate the curvature effect (this has to be con- 
firmed, however). 

Fifth, we have been talking about "typical" fluctuations, 
e.g., when we specified a one-cr amplitude of a given 
peculiar-velocity gradient invariant. Especially regional 
observational data reaching up to lOOMpc could well be 
those of an untypical region lying in the tail of the {non- 
Gaussian) probability distributions for the (second and 
third) invariants. 

In this line we would like to point out that, even if 
the fluctuations in number density (the first moment 
of the galaxy distribution) may not be present, fluctu- 
ations may show up in higher moments of the galaxy 
distribution. A recent investigation of subsets from the 
IRAS 1.2 Jy catalog revealed large fluctuations in the 
Minkowski functionals as well as the two-point correla- 
tions |i7| at least up to scales of 200ft~ 1 Mpc. The IRAS 
1.2 Jy catalog is an example of a sample that may look 
"homogeneous" as far as the number density of galaxies 
is concerned [Q, but features fluctuations already for 
the second moment and, most dramatic, for the higher 
moments. 

A number of open questions are to be considered. We 
already mentioned the quantification of the likelihood to 
find a region of a given size for given expectation val- 
ues of the invariants. The probability distributions of 
the averaged initial invariants would allow us to set like- 
lihood constraints on the cosmological parameters. An- 
other related task is to quantify cosmic variance: due to 
the scale-dependence of the "cosmological parameters" 
we can in principle infer their variations as a function 
of scale. Controlling these variations allows a statisti- 
cal assessment of how and if we approach a "scale of 
homogeneity" . Results from spatial variations of param- 
eters in an averaged inhomogeneous model may then be 
compared with observed fluctuations in galaxy or cluster 
catalogs. 



We have seen that the backreaction term, although 
it may be numerically small, significantly modifies the 
evolution of cosmological parameters, driven by the fluc- 
tuations in the velocity field. We expect exciting new 
insights by moving to general relativistic average mod- 
els, a subject of forthcoming work. 
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APPENDIX A: DERIVATION OF THE 
GENERALIZED FRIEDMANN EQUATION OF 
NEWTONIAN COSMOLOGY 

In this appendix we give a short derivation of the gen- 
eralized Friedmann equation [^ . 

In Eulerian space one of the dynamical fields is the 
velocity v(x, t). It is suitable to introduce the rate of 
expansion 6 = V ■ v, the shear (<7y ), and rotation uj = 
iVxv via a decomposition of the velocity gradient^: 

v i,j — a ij + 3%0 + u ij — ~ e i]kUk] = 0. 

(Al) 

The rate of shear a and the rate of rotation lu are defined 
as follows: 

a 2 = \ OijOij , uj 2 — \ LOijUJij. (A2) 

Now, consider the volume of an Eulerian spatial do- 
main T> at a given time, V = j v d 3 x, and follow the po- 
sition vectors x = f (X, t) of all fluid elements (labelled by 
the Lagrangian coordinates X within the domain). Then, 
each volume element changes according to d 3 x = Jd 3 A, 
where J is the determinant of the transformation from 
Eulerian to Lagrangian coordinates. The total rate of 
change of the volume of the same collection of fluid ele- 
ments of a domain may then be calculated as follows: 

= i / d 3 X 9J = (6) v = 3H V , (A3) 

where dt is the total (Lagrangian) time-derivative and 
Ht> = d-p/ar) the Hubble-parameter of the domain T>. 

It is crucial to notice that the total time-derivative 
does not commute with spatial averaging. For an arbi- 
trary tensor field A we derive with the help of the above 
definitions the following commutation rule: 

d t (A) v - (d t A) v = {A6) v - (6) v {A) v . (A4) 

The generalized Friedmann equation Eq. (Q) of the main 
text together with the backreaction term Eq. (||) follows 
by setting A = 6 and using the mass conservation Eq. (^) 
and Raychaudhuri's equation for the local evolution of 
the expansion rate 9: 

6 = A- 4ttG P --8 2 + 2(uj 2 - a 2 ). (A5) 
3 

(Eq. ( |A5| ) follows from the trace of the spatial derivative 
of Euler's equation d t v = g using the field equation V • 
g = A — AttGq for the gravitational field strength g.) 



6 A comma denotes partial derivative with respect to Eule- 
rian coordinates d/dxi =,i; 5ij is the Kronecker delta and 
tijk the totally antisymmetric unit tensor; we adopt the sum- 
mation convention. 
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APPENDIX B: INVARIANTS 

The three principal scalar invariants of a tensor field 
A = (Aij) in Cartesian coordinates are defined by 

= tr(Ay), 

II(A ij ) = ±(tr(A ij f-tr((A ij ) 2 )), 
III(^j) = det(^). 

In particular for A^ — Vij the invariants are expressible 
in terms of kinematical scalars and can be written as 
total divergences of vector fields, which has been used 
and discussed in the context of perturbation solutions 
MM: 



= Vi t i = V • v = 9, 
= |V- (v(V • v) - (v • V)v 



(Bl) 
(B2) 



, ,2 ^2 , l/)2 



= iV-(^V-(v(V-v)-(vV)v)v- 

(v(V • v) - (v • V)v) • Vv) (B3) 

= i0 3 + 20 ((T 2 + \US 2 ) + (JijCTjkCrki - CTijUJiLUj. 

These expressions are also valid for comoving quantities 
I(0g,Ui) = V 9 -u etc. , (B4) 
with the comoving differentials d qj = -£p and V ? . 

APPENDIX C: INITIAL CONDITIONS 

To estimate the magnitude of the volume-averages of 
the invariants, we first relate the initial peculiar-velocity 
field U(X) = u(q, ti) to the initial density contrast 
<5i(X) = (f?(X, ti) — qh (ti))/gn (ti) using the Eulerian lin- 
ear approximation. 

At time ti the comoving coordinates equal the La- 
grangian coordinates q = X by definition. Considering 
only the growing mode in an Einstein-de-Sitter back- 
ground <5(q, t) — (i/ii) 2 / 3 <5i(q), the peculiar-expansion 
rate evolves according to 8 u (q,t) = (t/tiY^O^q) with 



3 ti 



Si(d). 



(CI) 



Assuming a non-rotational initial velocity field 

k 



u(q,< ; ) 



d 3 A; 0i(h) — e ik q 

£r [ d 3 fc Ji(k)A e 



3 U 



ik 2 



(C2) 



with fc = |k| and the Fourier transform A(k) of a field 
•A(q): 



I(k) L d3q - 4(q)e ^ kq - 



(C3) 



In an Einstein-de-Sitter background the non-rotational 
initial velocity field is related to the initi al d isplacement 
field by ipuj = 3/2 iit/jU, and with Eq. ( |C2| ) we get the 
invariants of ip\ij in terms of the initial density contrast 
Si. 



1. Volume— averaged invariants 

To simplify calculations we consider spherical domains 
Br with radius R. Since we assume that the perturba- 
tions in the initial density field are stochastically inde- 
pendent from the specific location in space, we may cen- 
ter the sphere on the origin. We define the normalized 
window function (a top-hat) 



W R (X) 



i/\b r \ ifxetffl, 



o 



else, 



(C4) 



where \Br\ = 4ir/3 R 3 is the volume of the sphere. The 
Fourier transform of Wr, 



W R (k) 



\Br\ (2tt)2 k 3 



(sm(kR) - kRcos(kR)) , (C5) 



is depending on k = |k| only. The spatial average 
(^(X))^ of a field A may be written as a convolution: 

JR 3 

=(2tt) 3 / d 3 fc A(k)W R (k). (C6) 



With Eq. (|C2|), jBlD-jBal), and fl39|) we can express the 
averaged invariants in terms of the initial density contrast 
#i (remember U(X) = u(q,ii) and ip^j = 3/2 tiU^j): 



(C7) 



<I,) Sa = -(2tt) 3 / d 3 k Si(k)W R (k), 



d 3 fci / d 3 fc 2 



6i(ki)6i(k 2 ) Wfldkx +k 2 |) ( I - (k ^ 2 !g )2 ) (Cn 



h 2 h 2 



and 



(m,)» = -(2tt) 3 / d 3 ^ / d 3 k 2 [ d 3 k 3 

JR 3 JR 3 JR 3 

5(ki)5(k 2 )5(k3) ^fl(|k 1 + k 2 + k 3 |)F(k 1 ,k 2 ,k 3 ), 

(C9) 
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with 



F(ki,k 2 ,k3 



l(k2 
2 



h 2 h 2 



l(kj_ 

3 



k 2 )(k 1 -k 3 )(k 2 -k 3 ) 

k 2 k 2 k 2 



(CIO) 



Similar expressions without the spatial averaging were 



presented in a systematic way by 1 41 43] (see also 



To calculate the evolution of the scale factor using 
Eq. ( f4(| ) we also need the normalized over-density (S) v — 
((q)v ~ Qh)/qh- With Eq. (§ and g H (t) = g H (ti)/a(t) 3 
we can express {8) v in terms of the averaged first invari- 
ant: 



(1 + 



'V 



T>i) 



1 



(1 - (Ii) 



1. 



(Cll) 



2. Distribution of the volume— averaged invariants 



In Appendix C 1 we expressed the averaged invariants 
of the peculiar-velocity field as volume-averages over the 
initial density contrast Si. To get a handle on the distri- 
bution of these volume-averaged invariants we calculate 
their ensemble mean and their variance under the as- 
sumption that 5i is a Gaussian random field. A Gaus- 
sian density field is stochastically specified by its power- 
spectrum P\(k) with k = |k| (see e.g. JT^j), determining 
the correlation of the Fourier modes of the initial density 
contrast: 



E[tfi(k)tfi(k')] =S D (k + k') Pi(k). 



(C12) 



S D (-) is the Dirac delta-distribution. Contrary to the 
spatial average {')g a over a, finite domain, E[-] is denoting 
the ensemble average. Since a Gaussian field is ergodic 
p4j , a spatial average over the whole space {-}^3 equals 
the ensemble average E[-]. However, we focus on spatial 
averages over compact spherical domains Br. 

By the definition of a Gaussian random field all odd 
moments of the density contrast van ish. I.e. E[<5i(X) n ] = 
for n odd, and from Eqs. (|C7|-|C9|) follows 



E[(Ii> B J 
E[<IIIi}~ 



= = E[(I s ) Bfl 

= o = E[(ni)- 



:iii) B j> 



(C13) 



This already tells us that and (IIi) BR , as well as 

and are uncorrelated, and may be spec- 

ified independently. 

The other invariants and their products require some 
calculation s: 

With Eq. @, (|GL2|), and (|Cia| ) we obtain 



<x 2 (P)=E ((I i ) 6H -E[(I i ) BR ]) 2 = E[<I, 



(27r) b / d 3 /c Pi(k)W R {kY 



R 3 



(C14) 



The variance o"i(P) is equal to the well-known mean 
square fluctuations of the initial density contrast field 
smoothed with a top-hat of radius R. 
For the second invariant we obtain 



E[(IIi) B ] = 2(2tt) 3 / d 3 ^ Pi(fe) x 



d 6 k 2 5 u {k 2 ) 1 



(k!-k 2 -fc \f \~ 

fcf(k 2 - kl ^ r R{k2) - (C15) 



Using spherical coordinates and taking the limit k 2 
in the last integral, we arrive at 



E[(Hi) B J = 0. 



(C16) 



For a Gaussian density field the four-point correlation 
function factorizes into two-point correlations as ex- 
pressed by 

E[^ i (k 1 )J i (k 2 )5 i (k 3 )? i (k 4 )] = 

Pi(ki)Pi(h) 5 D (k 1 + k 2 )S D (k 3 + k 4 ) 
+P i (k 1 )P i (k 2 ) 5 D (k 1 +k 3 )5 D (k 2 +k 4 ) (C17) 
+P i (k 1 )P i (k 2 ) 5 D (kt + k 4 )S D (k 2 + k 3 ). 

Then, the varia nce of the second invariant can be calcu- 
lated from Eq. Q: 

a I 2 I (i?)=E[(II i ) e J = 



(27T) 



■ f d 3 fcx f d 3 k 2 Pi(fci)Pi(fc 2 ) 



WaQk! 



1 



(ki ■ k 2 

h 2 h 2 



A similar calculation gives 

E[(I,)„ (HI,)- ] = 0, 



(C18) 



(C19) 



showing that and are uncorrelated, and 

may be chosen independently for a Gaussian random 
field. 

Similar to Eq. ( [cil] ) one may factorize the six-point func- 
tion for a Gaussian density field. This enables us to ex- 
press E[(IIIi) 2 j H ] in terms of the power spectrum. After 
some lengthy algebra we obtain 



a I 2 n (i?)=E[(III i ) 2 J 



d 3 fc 2 



d 3 fe 



(27r) b / d 3 fci / u k 2 

Pi {h ) Pi (k 2 ) P ; (* 3 ) [w R {h ) 2 G x (^ , k 2 , k. 
M? fl (|k 1 + k 2 + k 3 |) 2 G 2 (k 1 ,k 2 ,k 3 ) 

with 

^(ki.ka.ka) = 
-2 



(C20) 
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-7> (C21) 
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and 



G a (ki,k 2l k3 



! (k 2 -k 3 ) 2 



2 (ki - kgHki -k 3 ) 

3 k-y k^ 

o (k 1 -k 2 ) 3 (k 1 -k 3 )(k 2 -k 3 ) 

1.4 1.4 1.2 
fi, 1 fi, 2 '<' 3 

2 (k 1 -k 2 )(k 1 -k 3 )(k 2 -k 3 ) 

3 h 2 k 2 h 2 

2(kl-k 2 ) 4 (k!-k 2 ) 2 

3" 



h 2 h 2 
K 1 K 2 



(C22) 



3. Some Numbers 

To calculate the fluctuations <Ji(R) etc. we need to 
specify the initial power spectrum Pi(k). Models for the 
linearly evolved power spectrum of the density contrast 
fluctuations at present time to are typically given by 



P(k,t )=Ak n T{k,t f 



(C23) 



with the amplitude A, the primordial spectral index 
n and the transfer function T(k, to). We use n = 1 
and a transfer function for Cold Dark Matter (CDM) 
parametrized according to Bardeen et al. p5[, for dark 
matter with il m = 1, a small baryon fraction, and three 
relativistic neutrino flavors: 

ln(l + rfc) 
T(Mo)= \k x 

(1 + sfc + (ifc) 2 + (ufc) 3 + (vkfy* , (C24) 

with r = 9.36Mpc, s = 15.56Mpc, t = 64.4Mpc, 
u = 21.84Mpc, and v = 26.84Mpc. Consistent with 
Sect. [V an Einstein-de-Sitter background with Ho 



50kms Mpc 1 is assumed, k is in units of Mpc 1 . We 
choose the normalization A = 2.19 10 4 Mpc 4 resulting in 
ct(8/0.5 Mpc) = 1 (see Eq. jC14| ) and Table |). This 
normalization is smaller than the COBE-normalization 
fhSf , but still larger than the cluster-normalization pTj] 
for this CDM-model. 

We choose our initial conditions at z\ = 200, since the 
linear power spectrum does not significantly change its 
shape for z < z\ [0] . The linear power spectrum at initial 
time t{, with a(ti) = 1, is dir ectly related to the linear 
power spectrum, given in Eq. ( C23| ) at present time to by 



Pi(k) = 



P(k,to) 
a(t Q ) 2 ' 



(C25) 



Using this CDM-model we calculate ai(R) and (Jn(R) 
for initial domains of radius R (see Fig. |l3| ). 
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FIG. 13. The values of ai(R) and an (R) are shown 
against the radius R of the spherical initial domain Br. On 
top also the radius a(to)R, with a(to) = 201, of the lin- 
early evolved domain is given. The dotted line is t he re- 
sult for oi(-R) using the truncated power of Eq. ( |C26[ ) . The 
short-dashed long-dashed line and the dashed-dotted line 
correspond to power spectra with a cut at small k, corre- 
sponding to L = 200Mpc and L = 400Mpc respectively. 

Some comments on the numerical methods are in or- 
der: calculating cr\(R) from Eq. (C14) does not impose 
any numerical problems. We used a crude Monte-Carlo. 



The triple integral appearing in Eq. (C18) for the calcula- 
tion of CTn(i?) required stratified sampling. In both cases 
our results are accurate to o ne di git at least. The cal- 
culation of oin(i?) from Eq. (C20) requires some major 
computational and algorithmic efforts in order to eval- 
uate the eight integrals. This is beyond the scope of 
this article. We only need an order of magnitude esti- 
mate for ctiii (R) . Guided by our analysis of the spheri- 
cal collapse model (see Eq. (p0[)) we will use the estimate 
cin(-R) ~ 1/27 <Ji(R) 3 in our calculations. The values 
for oi(i?), aii(-B), and am(R) used in our calculations 
in Sect. Ill are given in Table pi. 

It has been shown that the accuracy of the "Zel'dovich 
approximation" may be increased by using a smoothed 
initial density field, i.e. a truncated power spectrum P s 
g-|l|3]g]. The optimum smoothing scale for CDM 
initial conditions is in our case k s = 1.687MPC" 1 p2[ , 
resulting in 



Ps(k) 



-k 2 /k 2 B 



Pi(k). 



(C26) 



In Fig. [13| we show ai(R) calculated using P s (k) instead 
of Pi(k). Already for domains with a scaled radius larger 
than 20Mpc the difference becomes unimportant. 

In a periodic box with sidelength L, the power spec- 
trum is truncated at small k. The suppression of the 
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fluctuations in the initial conditions may be estimated 
by using a cut: P(k) = for k < 2n/L. In Fig. || the ef- 
fect of applying this cut to o\ (R) is illustrated for scaled 
sidelengths L = 200Mpc and L = 400Mpc, correspond- 
ing to spherical domains with a radius of aR — 125Mpc 
and aR = 250Mpc. 
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